close all 
clear all

rootdir = '/home/predatt/anatod/ISI';
datadir = '/home/predatt/anatod/ISI/raw'; 
gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG';
figdir = '/home/predatt/anatod/ISI/MEG/FIGURES/FINAL';

cd(gavgdir)
load neighbours
load subjname
load ChanSel
cond = 1:4;

clc

%% average over conditions of interest
load gavg_TFRs_TONES

cfg = [];
cfg.keepindividual = 'no';
cfg.foilim         = 'all';
cfg.toilim         = 'all';
cfg.channel        = 'MEG';
    
    strong_lo = ft_freqgrandaverage(cfg,gavg_lo{1}, gavg_lo{2});
    weak_lo = ft_freqgrandaverage(cfg,gavg_lo{3}, gavg_lo{4});
    att1_lo = ft_freqgrandaverage(cfg,gavg_lo{1}, gavg_lo{3});
    att2_lo = ft_freqgrandaverage(cfg,gavg_lo{2}, gavg_lo{4});

 % calculate differences for interactions

strong1_lo = gavg_lo{1};
strong2_lo = gavg_lo{2};
weak1_lo = gavg_lo{3};
weak2_lo = gavg_lo{4};

exp_diff_lo = gavg_lo{1};
unexp_diff_lo = gavg_lo{1};
att1_diff_lo = gavg_lo{1};
att2_diff_lo = gavg_lo{1};

exp_diff_lo.powspctrm = strong1_lo.powspctrm - strong2_lo.powspctrm;
unexp_diff_lo.powspctrm = weak1_lo.powspctrm - weak2_lo.powspctrm;
att1_diff_lo.powspctrm = strong1_lo.powspctrm - weak1_lo.powspctrm;
att2_diff_lo.powspctrm = strong2_lo.powspctrm - weak2_lo.powspctrm;
      
clc
%% Stats 

cfg = [];
cfg.avgoverchan = 'yes';
cfg.avgovertime = 'no';
cfg.avgoverfreq = 'no';
cfg.method = 'montecarlo';
cfg.statistic = 'depsamplesT';
cfg.numrandomization = 5000;
cfg.correctm = 'cluster';
cfg.neighbours = neighbours;
cfg.clusteralpha = 0.05;
cfg.clustertail = 0;
cfg.tail = 0;
cfg.alpha = 0.05;
cfg.correcttail = 'alpha';
Nsb = 1:length(subjname);
cfg.design = [ ones(1,length(Nsb)) 2*ones(1,length(Nsb)); ... 
    Nsb Nsb ]; 
cfg.ivar = 1; 
cfg.uvar = 2; 

% LOW frequencies
cfg.channel = channels_tfr_beta;
cfg.latency = [0 1.01];

    expectation = ft_freqstatistics(cfg,strong_lo,weak_lo); 
    attention = ft_freqstatistics(cfg,att1_lo,att2_lo);
    interaction = ft_freqstatistics(cfg,exp_diff_lo,unexp_diff_lo); 

clc
min(min(expectation.prob)) % 0.164
min(min(attention.prob)) % 0.290
min(min(interaction.prob)) % p=0.019

cfgoutput.alphalevel = 0.05;
tfr_output(cfgoutput,interaction)
 
% p=0.018, clusterstat=-269.45, freqs=13-32Hz, times=0.18-0.40sec 
% p=0.028, clusterstat=-200.80, freqs=15-35Hz, times=0.85-1.00sec 
% p=0.040, clusterstat=-170.43, freqs=15-29Hz, times=0.48-0.63sec 

figure(1);clf;imagesc(squeeze(interaction.negclusterslabelmat)==1)

clc
cfgoutput.alphalevel = 0.05;
tfr_output(cfgoutput, interaction);
% freqs=13-32Hz, times=0.18-0.40sec 

% post-hoc tests: constrain freq ant time to interaction cluster

cfg.avgovertime = 'yes';
cfg.avgoverfreq = 'yes';
cfg.latency = [0.1801 0.401];
cfg.frequency = [13 32];

    expectation_att1_cluster1 = ft_freqstatistics(cfg,weak1_lo,strong1_lo); 
    expectation_att2_cluster1 = ft_freqstatistics(cfg,weak2_lo,strong2_lo); 
    
cfg.latency = [0.85 1];
cfg.frequency = [15 35];

    expectation_att1_cluster2 = ft_freqstatistics(cfg,weak1_lo,strong1_lo); 
    expectation_att2_cluster2 = ft_freqstatistics(cfg,weak2_lo,strong2_lo); 

    cfg.latency = [0.48 0.631];
cfg.frequency = [15 29];

    expectation_att1_cluster3 = ft_freqstatistics(cfg,weak1_lo,strong1_lo); 
    expectation_att2_cluster3 = ft_freqstatistics(cfg,weak2_lo,strong2_lo);
    
   
clc
close all
%% display results

% p=0.018, clusterstat=-269.45, freqs=13-32Hz, times=0.18-0.40sec 
min(min(expectation_att1_cluster1.prob)) % 0.007, positive cluster
min(min(expectation_att2_cluster1.prob)) % non-significant

% p=0.028, clusterstat=-200.80, freqs=15-35Hz, times=0.85-1.00sec 
min(min(expectation_att1_cluster2.prob)) % nonsig
min(min(expectation_att2_cluster2.prob)) % nonsig

% p=0.040, clusterstat=-170.43, freqs=15-29Hz, times=0.48-0.63sec 
min(min(expectation_att1_cluster3.prob)) % 0.0046
min(min(expectation_att2_cluster3.prob)) % nonsig


close all
clc

%% regular ttests

cfgselect = [];
cfgselect.channel = channels_tfr_beta;
cfgselect.avgoverchan = 'yes';
cfgselect.avgoverfreq = 'yes';
cfgselect.avgovertime = 'yes';
cfgselect.latency = [0.18 0.401];
cfgselect.frequency = [13 32];

clust1_weak1 = ft_selectdata(cfgselect, weak1_lo);
clust1_weak2 = ft_selectdata(cfgselect, weak2_lo);
clust1_strong1 = ft_selectdata(cfgselect, strong1_lo);
clust1_strong2 = ft_selectdata(cfgselect, strong2_lo);

cfgselect.latency = [0.85 1.01];
cfgselect.frequency = [15 35];
clust2_weak1 = ft_selectdata(cfgselect, weak1_lo);
clust2_weak2 = ft_selectdata(cfgselect, weak2_lo);
clust2_strong1 = ft_selectdata(cfgselect, strong1_lo);
clust2_strong2 = ft_selectdata(cfgselect, strong2_lo);

cfgselect.latency = [0.48 0.631];
cfgselect.frequency = [15 29];
clust3_weak1 = ft_selectdata(cfgselect, weak1_lo);
clust3_weak2 = ft_selectdata(cfgselect, weak2_lo);
clust3_strong1 = ft_selectdata(cfgselect, strong1_lo);
clust3_strong2 = ft_selectdata(cfgselect, strong2_lo);

[h p ci stats] = ttest(clust1_weak1.powspctrm, clust1_strong1.powspctrm);
expectation_a1_cluster1.p = p;
expectation_a1_cluster1.tstat = stats.tstat;

[h p ci stats] = ttest(clust1_weak2.powspctrm, clust1_strong2.powspctrm);
expectation_a2_cluster1.p = p;
expectation_a2_cluster1.tstat = stats.tstat;

[h p ci stats] = ttest(clust2_weak1.powspctrm, clust2_strong1.powspctrm);
expectation_a1_cluster2.p = p;
expectation_a1_cluster2.tstat = stats.tstat;

[h p ci stats] = ttest(clust2_weak2.powspctrm, clust2_strong2.powspctrm);
expectation_a2_cluster2.p = p;
expectation_a2_cluster2.tstat = stats.tstat;

[h p ci stats] = ttest(clust3_weak1.powspctrm, clust3_strong1.powspctrm);
expectation_a1_cluster3.p = p;
expectation_a1_cluster3.tstat = stats.tstat;

[h p ci stats] = ttest(clust3_weak2.powspctrm, clust3_strong2.powspctrm);
expectation_a2_cluster3.p = p;
expectation_a2_cluster3.tstat = stats.tstat;
